      !=======================================================================
      !> @file aga-v1.f90
      !> @brief The implementation of AGA
      !> @authors Ary Rodriguez-Gonzalez & Alejandro Esquivel
      !> @date 23/Jun/2011
      !
      ! Copyright (c) 2011 Ary Rodriguez-Gonzalez & Alejandro Esquivel
      !
      ! This file is part of AGA-V1.
      !
      ! AGA-V1 is free software; you can redistribute it and/or modify
      ! it under the terms of the GNU General Public License as published by
      ! the Free Software Foundation; either version 3 of the License, or
      ! (at your option) any later version.
      !
      ! This program is distributed in the hope that it will be useful,
      ! but WITHOUT ANY WARRANTY; without even the implied warranty of
      ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
      ! GNU General Public License for more details.
      ! You should have received a copy of the GNU General Public License
      ! along with this program.  If not, see http://www.gnu.org/licenses/.
      !=======================================================================
      !> @brief Module containing the implementation of the algorithm
      !> @details It has to be called from a main program that serves as an 
      !! interface if used with a experimental/observational data.
      !> @todo Debug OpenMP implementation
      module agav1
      contains
      ! ----------------------------------------------------------------------
      !> @brief THE algorithm
      !> @details This is the implementation of the Asexual Genetic Algorithm
      !! (AGA), it requires that the experimental/observational data have been
      !! previously loaded (this is done in the main program, and was left out
      !! of the main algorithm to improve flexibility), and if MPI is enabled
      !! it also needs to be prevously loaded. The algorithm loads the rest of 
      !! the runtime parameters and return the fittest individual and an
      !! estimate of the errors.
      !> @param[in] xobs(ndata) : Vector of data abscissae
      !> @param[in] yobs(ndata) : Vector of data ordinates
      !> @param[in] ysigma(ndata) : Vector of uncertainties in the ordinate
      !> @param[out] q(npar) : Vector of parameters to be adjusted, at the end
      !! returns the fittest individul 
      !> @param[out] delta(npar) : Window (hyperbox) half-size, at the end
      !! the variable returns the uncertainty of q
      !> @param[in] ndata : Number of data points
        subroutine Algorithm(xobs,yobs,ysigma,q,delta,npar,ndata)
          use input                   ! located in extfunc.f90
          use func                    ! located in extfunc.f90
          use exfunc                  ! located in extfunc.f90
          use mpi_module              ! located in modules.f90
          use numgaus                 ! located in modules.f90
          implicit none
          integer :: npar, ndata, nsteps, nruns, nparent, nind, nerr,  &         
                     deltaconst, nmig, ngaussdev, igaus, gensize
          real :: beta,reduc,fito,deltafit
          real :: mean, sigma
          !
          real, dimension(:), allocatable :: q0,delta0, q, delta, fit
          real, dimension(:), allocatable :: xobs,yobs,ysigma
          !
          real, dimension(:,:),allocatable :: arrmain,arrpar,arrgen,arrerr,  &   
                                              qbounds
          !
          integer :: i,j,k,niter,idum,nit,ij,ncrit
          integer :: stk, nproc, omp_get_thread_num, KMP_GET_STACKSIZE_S,nprocs
          integer :: omp_get_num_threads
          !
          integer :: rand_size
          integer, allocatable, dimension(:) :: rand_seed
          character (len=10) :: system_time
          real :: rtime
          !
#ifdef MPIP
          real, dimension(:,:),allocatable :: arrerrtot
#endif
          !
          !     OpenMP initialization
          !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
          !$    WRITE (6,*) 'Compiled with OpenMP'
          !$OMP PARALLEL PRIVATE(NPROC)
          !$           stk = KMP_GET_STACKSIZE_S()
          !$           nproc = OMP_GET_THREAD_NUM()
          !$           WRITE(6,*) 'Processor (thread) #:',nproc,' ready'
          !$OMP END PARALLEL
          !$           WRITE(6,*) 'Stack size per thread : ', stk/1048576,' Mb'
          !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
          !
          !   initializes seed for random number generator
          
          call random_seed(size=rand_size)
          allocate(rand_seed(1:rand_size))
          call date_and_time(time=system_time)
          read(system_time,*) rtime
          rand_seed=int(rtime*1000.)
#ifdef MPIP
          rand_seed=rand_seed*rank
#endif
          call random_seed(put=rand_seed)    
          deallocate(rand_seed)
          !----------------------------------------------
          !
#ifdef MPIP
          do i=0,np-1
                      if (rank.eq.i) then
#endif
                              call InitCond(nsteps,nruns,nparent,nind,nerr, deltaconst,   &
                                                    nmig, ngaussdev, igaus, npar, gensize, reduc, &
                                                beta, qbounds,q0,delta0,fit,arrmain,arrpar,arrgen )                     
              
                              allocate(arrerr(npar,nerr/np),q(npar),delta(npar))
                              q(:)=q0
                              delta(:)=delta0
#ifdef MPIP
                              allocate(arrerrtot(npar,nerr))
             endif
             call mpi_barrier (mpi_comm_world,err)
          enddo
#endif
          !
          call spawn(npar,nparent,nmig,ngaussdev,delta,q,qbounds,arrpar)
          !   
          ! Main loop
          !  
          do j=1,nerr/np
                Print *,'Nerr=',j,np
             !
              if (nerr .ne. 1) call randdata(ndata,yobs,ysigma)
             !
             nit=1
             !
             do while(nit.le.nruns)
                !
                deltafit=1e20
                fit(1)=1e20
                niter = 1
                ncrit=0
                !
                do while(niter .le. nsteps)
       
                  call deltabox(npar,niter,nind,deltaconst,reduc,arrmain,     &
                                delta0,delta)         
                   !
                   !$omp parallel shared(arrgen,arrmain,arrpar,fit) &
                   !$omp           private(nproc,nprocs)
                   !$ nproc = OMP_GET_THREAD_NUM()
                   !$ nprocs= omp_get_num_threads()
                   !$omp do SCHEDULE(static,nparent/nprocs)
                   !
                   do i=1,nparent
                      call spawn(npar,gensize,nmig,ngaussdev,delta,           &
                                 arrpar(:,i), qbounds, arrgen )
                      !
                      arrmain(:,((i-1)*gensize+1):i*gensize)=arrgen(:,:)
                   enddo
                   !$omp end do
                   !
                   !$omp do SCHEDULE(static,nind/nprocs)
                   do i=1,nind
                      !!   this works only for sums of gaussians
                      If (igaus .eq. 1) Then
                         !
                         call gaussian(ngaus,i,arrmain)
                         !
                      EndIf
                      q(:)=arrmain(:,i)
                      fit(i)=meritfunc(xobs,yobs,ysigma,arrmain(:,i),         &
                                       npar,ndata)
                   enddo
                   !$omp end do
                   !$omp end parallel
                   !
                   call pickbest(npar,nind,nparent,arrmain,fit,arrpar)
      
                   niter = niter+1
                   deltafit=abs(1-fito/fit(1))
                   fito=fit(1)
      ! CRITERIO NUEVO
      !             call criterion(nparn,q,delta,beta,ncrit)
      !             print*,nit,niter,ncrit,j
      !             if (ncrit .eq. 1) exit
      !
      !             if (deltafit .lt. beta) exit
                    print*,niter,nit
                    if (fit(1) .lt. beta) exit
      !
                enddo
                nit=nit+1
                !
             enddo
            arrerr(:,j)=arrpar(:,1) !c(:)
          end do
          !
#ifdef MPIP
          call mpi_gather (arrerr,    npar*nerr/np, mpi_real8 ,               &
                           arrerrtot, npar*nerr/np, mpi_real8 ,               &
                            master, mpi_comm_world, err )
          if (rank.eq.master) then
#endif
          !
          open(2,file='coeffit.out')
          !
          do i=1,npar
#ifndef MPIP
             q(i)=med(arrerr(i,:),nerr)
             delta(i)=sigm(arrerr(i,:),nerr)
#else
             q(i)=med(arrerrtot(i,:),nerr)
             delta(i)=sigm(arrerrtot(i,:),nerr)
#endif
          enddo
          if (rank.eq.master) then
             do i=1,npar
                write(2,*)q(i),delta(i)
             enddo
                write(2,*)fit(1)
          endif
          close(2)
#ifdef MPIP
       endif
       dellocate(arrerrtot)
       call mpi_barrier (mpi_comm_world, err)
#endif
       deallocate(fit,arrerr)
       
      end subroutine Algorithm
      end module agav1
      !-----------------------------------------------------------------------
